
#include "fluid.h"
#include <iomanip>
#include <iostream>
using namespace std;

Fluid::Fluid(const Config &cfg, int imax, int jmax)
    : equsType{cfg.GetEquationType()}, isRestart{cfg.GetRestartTag()},
      nondimType{cfg.GetNondimType()}, reNum{cfg.GetReynoldsNumber()},
      mach{cfg.GetMachNumber()}, aoa{cfg.GetAngleAttack() * RAD},
      Tinf_d{cfg.GetTemperatureInf()}, reLen{cfg.GetReynoldsLength()},
      pInf_d{cfg.GetPressureInf()}, imax{imax}, jmax{jmax}

{
  int imaxp1{imax + 1}, jmaxp1{jmax + 1};
  cv.ReConstruct(imaxp1, jmaxp1);
  press.ReConstruct(imaxp1, jmaxp1);

  velInf_d = mach * sqrt(gama * rgas_d * Tinf_d);
  uInf_d = velInf_d * cos(aoa);
  vInf_d = velInf_d * sin(aoa);

  if (equsType == Equation_t::Euler)
  {
    dvE.ReConstruct(imaxp1, jmaxp1);
    rhoInf_d = pInf_d / (rgas_d * Tinf_d);
  }
  else
  {
    dvNS.ReConstruct(imaxp1, jmaxp1);
    muInf_d = SmuRef_d * pow(Tinf_d / STref_d, 1.5) * (STref_d + Sconst_d) /
              (Tinf_d + Sconst_d);

    rhoInf_d = reNum * muInf_d / (velInf_d * reLen);
    pInf_d = rhoInf_d * rgas_d * Tinf_d;
  }

  // 后面所有的物理量均是基于下面的变量 统一接口
  uInf = uInf_d;
  vInf = vInf_d;
  velInf = velInf_d;
  rhoInf = rhoInf_d;
  pInf = pInf_d;
  Tinf = Tinf_d;
  cpgas = cpgas_d;
  rgas = rgas_d;
  dynamicPress = 0.5 * rhoInf * velInf * velInf;

  sutMu = SmuRef_d;
  sutT0 = STref_d;
  sutTs = Sconst_d;
  muInf = muInf_d;
}

void Fluid::PrintEuler() const
{
  // cout << defaultfloat;
  cout << "mach = " << mach << "\n"
       << "temperature = " << Tinf_d << "\n"
       << "free stream rho = " << rhoInf_d << "\n"
       << "free stream velocity = " << velInf_d << "\n"
       << "free stream pressure = " << pInf_d << "\n"
       << "angle of atanck = " << aoa / RAD << "  [degree]\n"
       << "cpgas = " << cpgas_d << "\n"
       << "cvgas = " << cvgas_d << "\n"
       << "gas constant " << rgas_d << "\n";
}

void Fluid::PrintNS(void) const
{
  cout << "mach = " << mach << "\n"
       << "Reynolds = " << reNum << "\n"
       << "free stream temperature = " << Tinf_d << "\n"
       << "free stream rho = " << rhoInf_d << "\n"
       << "free stream velocity = " << velInf_d << "\n"
       << "free stream pressure = " << pInf_d << "\n"
       << "free stream mu = " << muInf_d << "\n"
       << "angle of atanck = " << aoa / RAD << "  [degree]\n"
       << "reynolds length = " << reLen << "\n"
       << "cpgas = " << cpgas_d << "\n"
       << "cvgas = " << cvgas_d << "\n"
       << "gas constant " << rgas_d << "\n"
       << "Sutherland:\n"
       << "muref = " << SmuRef_d << "\n"
       << "Tref = " << STref_d << "\n"
       << "const = " << Sconst_d << "\n";
}

// 暂时不初始化依赖变量
void Fluid::InitFlow()
{
  ConsVar_t initConv{rhoInf, rhoInf * uInf, rhoInf * vInf, 0.0};
  initConv.rhoE = pInf / gam1 + 0.5 * rhoInf * velInf * velInf;
  cv = initConv;
}

void Fluid::DependentVarsAllEuler(void)
{
  real_t pij;
  for (int j = 0; j <= jmax; ++j)
    for (int i = 0; i <= imax; ++i)
    {
      const ConsVar_t &conv{cv(i, j)};
      DeptVarE_t &dpv{dvE(i, j)};
      press(i, j) = pij = gam1 * (conv.rhoE -
                                  0.5 * (conv.rhou * conv.rhou + conv.rhov * conv.rhov) /
                                      conv.rho);
      dpv.T = pij / (rgas * conv.rho);
      dpv.c = sqrt(gama * rgas * dpv.T);
    }
}

void Fluid::DependVarsAllNS(void)
{
  for (int j = 0; j <= jmax; ++j)
    for (int i = 0; i <= imax; ++i)
    {
      const ConsVar_t &conv{cv(i, j)};
      DeptVarNS_t &dpv{dvNS(i, j)};
      press(i, j) = gam1 * (conv.rhoE -
                            0.5 * (conv.rhou * conv.rhou + conv.rhov * conv.rhov) /
                                conv.rho);
      dpv.T = press(i, j) / (rgas * conv.rho);
      dpv.c = sqrt(gama * rgas * dpv.T);
      dpv.mu = sutMu * pow(dpv.T / sutT0, 1.5) *
               (sutT0 + sutTs) / (dpv.T + sutTs);
      dpv.ka = dpv.mu * cpgas / Prl;
    }
}